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ABSTRACT 

We analyze and compare the bulges of a sample of L» spiral galaxies in hydro dynamical 
simulations in a cosmological context, using two different codes, P-DEVA and GASOLINE. The 
codes regulate star formation in very different ways, with P-DEVA simulations inputing low star 
formation efficiency under the assumption that feedback occurs on subgrid scales, while the 
GASOLINE simulations have feedback which drives large scale outflows. In all cases, the marked 
kncc-shapc in mass aggregation tracks, corresponding to the transition from an early phase of 
rapid mass assembly to a later slower one, separates the properties of two populations within 
the simulated bulges. The bulges analyzed show an important early starburst resulting from the 
collapse-like fast phase of mass assembly, followed by a second phase with lower star formation, 
driven by a variety of processes such as disk instabilities and/or mergers. Classifying bulge stellar 
particles identified at ^ = into old and young according to these two phases, we found bulge 
stellar sub-populations with distinct kinematics, shapes, stellar ages and metal contents. The 
young components are more oblate, generally smaller, more rotationally supported, with higher 
metallicity and less alpha-element enhanced than the old ones. These results are consistent 
with the current observational status of bulges, and provide an explanation for some apparently 
paradoxical observations, such as bulge rejuvenation and metal-content gradients observed. Our 
results suggest that bulges of galaxies will generically have two bulge populations which can 
be likened to classical and pseudo-bulges, with differences being in the relative proportions of the 
two, which may vary due to galaxy mass and specific mass accretion and merger histories. 

Subject headings: galaxies: bulges, galaxies: formation, galaxies: star formation, galax;ies: kinematics 
and dynamics, cosmology: theory, methods: numerical 
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1. Introduction 

Understanding how bulges form and evolve 
is important within the theories and models of 
galaxy formation and evolution. Bulges are the 
component responsible for the central light in ex- 



cess of the exponential disk (Freeman 1970), and 
account for more than 25% of the starlight emit- 
ted in the local Universe. There are three main 
classes of observed bulge: classical, pseudo-, and 



boxy or peanut-shaped bulges, see Kormendy fc 



Kennicutt (2004) and Athanassoula (2005) for a 
detailed discussion on bulge classification. The 
distinction is important in that different classes 
cloud have formed and evolved in different ways. 

The Milky Way Bulge is the only case in which 
individual stars are resolved, and thus provides 
unique information about bulge properties. The 
Milky Way is considered to have a boxy bulge, yet 
increased evidence for an old, a-enriched stellar 
population formed on a short time-scale has re- 



sulted in a two-component model (e.g. Tsujimoto 
& Bekki 2012) of the Bulge. It had been shown 



that two stellar populations coexist in the Bulge, 
separated in age and metallicity (McWilliam & 



Rich 


1994| Feltzing & Gilmore|2000, ,Barbuy et al. 


1999 


van Loon et al.|2003| 


Groenewegen & Blom- 


maert 2005 Zoccali et al. 


2006 Fulbright et al. 


2007 ZoccaU et al. 2008) 


and that the separa- 



tion somewhat extends to kinematics ( Zhao et al. 



1994 Soto et al. 2007), even if age determina- 



tions through color- magnitude diagrams showed 
that most bulge stars in the Galaxy are older than 
10 Gyrs (Ortolani et al. 1995; Feltzing & Gilmore 
2000} [Zoccah et al.||2006| ,Clarkson et al.||2008[ ). 



Recently Bensby & et al. (20111 reported a 



2-13 Gyrs age span among microlensed turn-off 
Bulge stars. These latest data from high resolu- 
tion spectrography, have put increased emphasis 
on the picture of a two-component Bulge. Indeed, 
analyses of giant stars confirmed the presence of 
two distinct populations: a metal-poor enriched 
in [a/Fe] one with kinematics consistent with an 
old spheroid, and a metal-rich one with roughly so- 



lar [a/Fe] and with bar-like kinematics (Babusiaux 
et al.||2010| |Hill et al.|[20TT] [Gonzalez et al.||2011 



De Prop ris et al.[2011 Johnson et al.||2011 Robin 
et al. 2012nSoto et al.|2012||Ness fc Freemaii|20T2 



Ness et al. 2012). Comparing metallicities and 
compositions at different galactic latitudes, it has 



been found that while the [a/Fe] remains roughly 
constant and metal-poor stars show a remarkable 



homogeneity in 


the Bulge (|Lecureur et al.| 


2007 


Johnson et al. 


2011 




Gonzalez et al. 


2011 


, the 



most metal-rich stars near the galactic plane dis- 
appear at higher latitudes (see Zoccali et al.|2008" 



[Lecureur et al.|2007| [Babusiaux et al.|2010| among 
others). Further, the younger population is asso- 
ciated with the bar (Babusiaux et al. 2010 and 



references therein). These differences could be an 
indication that younger, more metal-rich stars in 
the Bulge define a smaller region than metal-poor 



ones (see also Robin et al. 2012 ). Otherwise, John- 



son et al. (2012) explored the link between Bulge 
and thick disk formation and found diverging be- 
haviors in the cases of [Na/Fe] and [La/Fe], which 
results in the Bulge resembling more the spheroid. 

In the bulges of external galaxies, similar con- 
clusions have been reached. Stellar population 
studies suggest that in many cases a secondary 
stellar population is superimposed on an older one 



(Ellis et al. 2001 Thomas fc Davies 2006 Car- 



oUo et al. 2007), and that the two populations 
are kinematically distinguishable, with the old 
population having spheroid-like kinematics, while 
this secondary population is more disk-like, see 



for example Prugniel et al. (2001), and the re- 



sults from the SAURON project (Peletier et al 



2007 Erwin 2008). Bulge stellar masses are gen- 
erally dominated by the old populations, with 
the young ones contributing less than a 25% in 
most cases ( MacArthur et al.||2009 ), although 'Ko-| 



rmendy et al. (2010) points out that significant 
numbers of local massive spiral galaxies appear to 
have dominant pseudo-bulges. 



Jablonka et al. (2007) conclude that most ex- 



ternal bulges in their observed sample are more 
metal rich and have lower [a/Fe] enhancement in 
their central regions than in their outer parts, a re- 



sult consistent with Moorthy & Holtzman ( 2006 1 , 
and with Milky Way Bulge results. MacArthur] 



et al. (2009) find a wide range of gradients, both 
positive and negative, allowing for different bulge 
formation mechanisms. A complementary piece of 
information are age gradients, where most authors 
find that the central regions of external bulges 
are younger than the outer ones, (e.g. 
Hortzman„2006{ [Jablonka et al.[|2007 



et al. 2009 Sanchez-Blazquez et al. 2011). Note, 



Moorthy & 



MacArthur 



however, that the last authors have also found neg- 
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ative age gradients in some cases. 

Summarizing, it would appear that bulge for- 
mation and evolution in a cosmological context 
has to account for a duality in stellar popula- 
tions, and for the coexistence of spheroid-like fea- 
tures with disk-like ones, whilst some particular- 
ities may be related to minor merger events (see 
ussion in Combes 20091. Theoretical models fo- 
cusing on the physical processes responsible for 
bulge properties have a long history, a sample of 
which is presented below. Generally, the classical 
type is thought to originate from fast gas collapse 
at high redshift, or from gas clumps in a proto-disk 
that are driven to the central regions by dynamical 
friction. Whilst instabilities of the disk, perhaps 
associated with bars, are believed to result in the 
formation of pseudo- and/or boxy bulges. 

Metal enrichment in bulges was first analyt- 
ically studied through pure chemical evolution 
models. Matteucci & Brocato (1990) first pre- 



dicted an [a/Fe] enhancement at bulges as a con- 
sequence of the different SNe I and II nucleosyn- 
thetic time-scales, observationally confirmed by 
McWilliam & Richl ^199^. IMolla et al.l (I2OOOI) 



used a multiphase multi-zone chemical evolution 
model and provided solutions to bulge chemical 
abundances and spectral indices. Ferreras et al. 



( 2003 ) found that very short infall time-scales are 
required for bulges. These models provided the 
basis of later developments to be implemented 
in hydrodynamical simulations, and, at the same 
time, gave us clear insights into some of the phys- 
ical processes involved. 

On the other hand, purely collisionless numer- 
ical studies of bulge formation have focused on 
morphology and dynamics by analyzing disk evo- 
lution in pre-prepared simulations. Such stud- 
ies showed that spheroidal-like bulges cannot be 



formed through bar-buckling instabilities (Debat 
tista et a"L||2004 1 while disk-like and boxy/peanut 



bulges are different in their nature, although both 
classes are associated with the presence of a bar 



( Athanassoula 


2005 


Moral et al. 


( 


2011) 



2005[ ). On the other hand, [Eliche- 
2011) analyze the effect of minor 
mergers on the inner part of disk galaxies, find- 
ing this process to be efficient in forming rotation- 
ally supported stellar inner components, i.e. disks, 



rings or spiral patterns. Hopkins et al. (20101 use 



a different approach, by constructing semi empir- 
ical models, based on observationally motivated 



halo occupation numbers from the Millenium sim- 
ulation, aiming to quantify the relative effect of 
galaxy mergers on bulge formation. They find ma- 
jor mergers to be the dominant mechanism for 
bulge and spheroid formation and assembly, while 
minor mergers contribute relatively more in lower 
mass systems. 

The first numerical studies of bulge formation 
using dissipative collapse by |Samland fc Gerhard] 
(2003) model the formation of a large disk inside 
a spinning (A — 0.05), growing dark matter (DM) 
halo, with an added accretion history taken from 
the large-scale simulations of the GIF- VIRGO 
consortium (Kauffmann et al. 1999). Chemical 



evolution is followed through two fiducial elements 
tracing the fraction of heavy elements produced 
by type la or type II SNe. The resulting bulge 
consists of at least two stellar sub-populations, an 
early collapse population and another one that 
formed later in the bar. INakasato fc Nomotol 



(2003) presented results of the evolution of a 
spherical Bcr top-hat over-dense region of 1.4 Mpc 
co-moving radius in rigid rotation evolved with a 
gravohydrodynamical code. Again the nucleosyn- 
thetic yields of type la or type II SNe are used, 
tracing the chemical enrichment in metallicity, O 
and Fe. Their results suggest that bulges consist 
of two chemically different components; one that 
has formed quickly through a sub-galactic merger 
in the proto-galaxy, and another one that formed 
gradually in the inner disk. |Kobayashi fc Nakasato] 
(2011) further developed this latter chemical evo- 



lution implementation, and used it to simulate 
a Milky Way-like galaxy. Their kinematic and 
chemical results follow closely the observed prop- 
erties of the Galaxy halo, bulge and thick disk. 
These simulations, however, do not take into ac- 
count the cosmological gas infall. 

Few studies on bulge formation have been made 
within a fully cosmological context. [Tissera fc 



Dominguez-Tenreiro ( 1998 1 and Governato et al 



( 2009 ) have studied the effects of mergers on classi- 



cal bulge stellar populations. Guedes et al. (2011) 
have run the highest resolution up-to-now simula- 
tion of Milky Way like disk formation, and got re- 
alistic disk and bulge properties but their focus on 
the details of bulge formation is relatively minor. 



Okamoto (2012) finds that pseudo bulges form, in 
his simulated Milky Way galaxies, by rapid gas 
supply at high-redshift, with their progenitors ob- 
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servable as high-redshift disks, and that this oc- 
curs prior to formation of the final disk. |Brook| 



et al. (2012a I obtained a lower mass late-type 
disk galaxy which has a bulge that grows from 
z = 1 mainly through purely secular processes. 
Other authors have put more emphasis on metal 
enrichment. For example, iRahimi et al. (2010) 



For example 
run a fully cosmological simulation with the GCD+ 
code ( Kawata fc Gibson|2003 ) , which incorporates 



chemical enrichment both by SNe la (Iwamoto 



et al. 


1999 


Kobayashi et al. 


2000 


(Woosley & Weaver 


1995 


I, as well 



and SNe II 



from intermediate mass stars. The code does not 
include a mechanism to diffuse metals between gas 
particles, resulting in an artificially high spread in 
the metallicity distributions, but robust averages. 
Their results underline the importance of mergers 
in bulge formation and their possible kinematic 
implications, the dependence of metal content on 
age, and the existence of accreted stars within the 
bulge. 

Thus it would appear that a variety of conclu- 
sions are being drawn from different groups using 
different codes and different physical models for 
their galaxy formation simulations. By analyz- 
ing in a unified manner (measuring the relevant 
properties with the same pipeline) simulated disk 
galaxies that are run with different codes and dif- 
ferent physical prescriptions and merger histories, 
we hope in this study to shed light on what pro- 
cesses of bulge formation, and subsequent bulge 
properties, are common within such simulations. 
Specifically, we analyze the bulges of three of the 



Moral et al. 


(20121, run with the P-DEVA (Serna 


et al. 


2003 


Martinez-Serrano et al. 


2008|) code, 



and those of two GASOLINE galaxies described in 
Brook et ai] ( |2012a[ ) and G. Stinson et al. (2012, 
submitted), with the aim of deciphering the pat- 
terns of bulge formation by focusing in the prop- 
erties of their stellar populations. 



Domenech-Moral et al. (2012) analyze disk for- 



mation in a cosmological context by running zoom 



simulations with P-DEVA. This code (Martinez 



Serrano et al. 2008 1 incorporates a statistical im- 
plementation of chemical enrichment based on |Tal- 
bot & Arnett (19731, including chemical feedback 



both by SNe la and SNe II, as well as mass loss 
from intermediate mass stars, involving 11 ele- 
ments. Radiative cooling takes into consideration 



the full element distribution at each point and 
time through a particular implementation of the 
Dimension Reduction regression ( |Weisberg|[2002| , 
while a SPH metal diffusion term mimics turbulent 
effects (Monaghan [2005 ). They have produced 
disk systems whose different components (thin and 
thick disk, halo, bulge) have properties nicely con- 
sistent with observations, for example, they have 
g-band B/T ratios between 0.13 and 0.36. In par- 
ticular, bulge metallicity and [a/Fe] distributions 
show bimodal patterns, that they interpret as re- 
sulting from fast and slow modes of star formation. 
The simulations run with P-DEVA assume that su- 
pernova feedback works on sub-grid scales and re- 
sults in a low star formation efficiency, which is 
thus used as an input parameter (see discussion 
in Agertz et al. 2011) that implicitly mimics the 
energetic feedback. 

On the other hand, the simulations run with the 
GASOLINE code have explicit feedback from mas- 
sive stars which drives large scale outflows. The 
two particular GASOLINE simulations were chosen 
for this study because (i) G-1578411 was shown 
to have a bulge which is mostly formed after the 
merger epoch. It is the late type simulated disk 



galaxy from Brook et al. (2012a I. A secular bulge 
grows between z = 1 and 2 = 0, driven at least 
partly by a bar. The final bulge to total light ra- 
tio is B/T=0.21. (ii) G-1536 is the simulated L* 
galaxy run by our MaGICC program that matches 
the widest range of observed galaxy properties (see 
G. Stinson et al. (20 12, submitted) and SG5LR in 
Brook et al.] ( |2012b[ )). It has a B/T ratio of 0.35. 



We use this suite of 5 inhomogeneous simulated 
disk galaxies to analyze the mass-weighted three- 
dimensional shape, size and kinematics of the 
bulges, as well as their mass-weighted age, metal- 
licity and chemical composition. More specifically, 
we first show that the simulated bulges consist of 
roughly two stellar populations (old and young) 
whose properties are correlated with their shape, 
kinematics, metallicity and composition, in line 
with the observational results we have previously 
described. In view of this agreement with observa- 
tions, we investigate the physical processes under- 
lying bulge formation and particularly, their rela- 
tion with the dynamics of the cosmic-web at high 
rcdshift, and of their host galaxies at low zs. 

An interesting aspect of this work is to clarify to 
what extent the formation mechanisms proposed 
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for massive ellipticals also apply to the formation 
of bulges, as suggested by similarities in observed 



properties (e.g. 


Franx 


1993 


Peletier et al. 


1999 


Carollo et al. 2001 


Bureau 2002 


Peletier 20081. In 



this regard, let us recall that the characteristics of 
mass assembly and star formation rate histories in 
ellipticals can be interpreted in terms of dark halo 
dynamics and its consequences. Indeed, analytical 
models as well as N-body simulations show that 
two different phases can be distinguished along 



the halo mass assembly process (Wechsler et al 



2002 Zhao et al. 2003): i) first a violent fast 



phase with high mass aggregation rates, resulting 
from collapse-like and merger events, and ii) later 
on a slow phase with much lower mass aggrega- 
tion rates. Small box hydrodynamical simulations 



( Dommguez-Tenreiro et al.|2006 


I as well as larger 


box ones ( 


Oser et al. 2010 


Dommguez-Tenreiro 


et al. 


2011 


1 confirmed this scenario as well as its 



implications on elliptical properties at low z { Cook 



et al. 2009). This scenario nicely explains appar- 
ently paradoxical observational data of elliptical 
galaxies. 

The paper is structured as follows. In Section[2] 
we describe the details of the simulations. The 
criteria used to select bulge stars are discussed in 
Section[3l In Section|4] we construct the star for- 
mation histories of the bulges and define the selec- 
tion criteria for the old and young stellar compo- 
nents. In the Scction[5]we discuss their shapes and 
kinematics, and in Section[6]thcir chemical compo- 
sition. The origin of the bulge stars in the frame- 
work of host galaxy and cosmic-web dynamics is 
given in Sectionj?) Finally, in SectionjS] we draw 
our conclusions after a brief summary of our find- 
ings and discuss the bulge formation scenario that 
emerges from our simulations. 

2. The simulations 

The main properties and resolution of our five 
simulated galaxies are shown in Table [T] All 
galaxies have previously appeared in the litera- 
ture, where more details can be garnered. Here 
we outline the main features of the codes and sim- 
ulations. 

2.1. P-DEVA 

We use the OpenMP parallel version of the 
DEVA code ( Serna et al.||2003 ), which includes the 



chemical feedback and cooling methods described 
in Martinez-Serrano et al. ( |2008 ), and in which 
the conservation laws (e.g. momentum, energy, 
angular momentum and entropy) hold accurately 



(Serna et al. 


20031. Full details of the 


recipes 


used are found in 


Domenech-Moral et al. 


( 


2012) 



where the simulations are first presented. The star 
formation recipe follows a Kennincutt - Schmidt- 
like law with a given density threshold, p*, and 
star formation efficiency . In line with |Agertz| 
et al. (2011), we implement inefficient SF param- 
eters (see Table [ij , which implicitly account for 
the regulation of star formation by feedback en- 
ergy processes by mimicking their effects, which 
are assumed to work on sub-grid scales. 



The chemical evolution implementation (Martinez 



Serrano et al. 2008) accounts for the dependence 
of radiative cooling on the detailed metal composi- 
tion of the gas, by means of a fast algorithm based 
on a metallicity parameter, C{T)- The code also 
tracks the full dependence of metal production 
on the detailed chemical composition of stellar 



particles (Talbot & Arnett 1973), through a Q 



formalism implementation of the stellar yields, 
for the first time in SPH. The delayed gas resti- 
tution from stars has been implemented through 
a probabilistic approach, that reduces statistical 
noise when compared with previous approaches, 
and therefore allows for a fair description of ele- 
ment enrichment at a lower computational cost. 
Moreover, the metals are diffused in such a way as 
to mimic the turbulent mixing in the interstellar 
medium. 

The simulations use the cosmological "zoom- 
in" technique, with high-resolution gas and dark 
matter in the region of the main object. The cos- 
mological parameters of a ACDM model were as- 
sumed (Qa = 0.723, n,n = 0.277, ilt = 0.04, and 
h = 0.7), in a 10 Mpc per side periodic box. Stellar 
masses are distributed according to the Salpeter 
initial mass function (IMF) ( |Salpeter|[T955l ), with 
a mass range of [M;,M„] = [0.1,100]Mq. 

2.2. GASOLINE 



The GASOLINE galaxies are cosmological zoom 
simulations derived from the McMaster Unbiased 
Galaxy Simulations (MUGS, [Stlnson et al.||2010 |. 
In G-1578411, the initial conditions (ICs) are 
scaled down, so that rather than residing in a 68 
h~^ Mpc cube, it is inside a cube with 34 h^^ 
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Object 


SMbar 




p* 




IMF 




^bulge 




Mgas 






IO^Mq 


h^-'^kpc 


cm^'^ 


% 




10^^ ergs 


kpc 


lOi°M0 


lO^M© 


Mpc 


G-1578411 


0.2 


0.15 


9.4 


1.7 


Kroupa93 


1.0 


1.30 


0.83 


1.64 


34 


G-1536 


1.9 


0.15 


9.4 


3.3 


Chab03 


1.0 


2.10 


0.96 


5.92 


68 


HD-5004A 


3.94 


0.2 


6 


1.0 


Salp55 




1.00 


1.50 


3.67 


10 


HD-5004B 


3.94 


0.2 


10 


0.8 


Salp55 




1.55 


1.60 


5.33 


10 


HD-5103B 


3.78 


0.2 


12 


0.8 


Salp55 




1.73 


1.47 


4.61 


10 



Table 1: The initial mass of gas particles {SMbar), minimum SPH smoothing length (hjo/t), density 
threshold (p*), star formation efficiency (c*), initial mass function (IMF), SN feedback (Esjv)j bulge 
radius (ibuige), stellar (M*) and gas (Mgas) mass of the simulated bulges (at z — 0), and periodic box 
length (Lbox)- The mass values correspond to the position selection. The P-DEVA runs use a fixed mass 
for the baryonic particles, equal to SMbar- 



Mpc sides, while G-1536 uses the same ICs as in 
the MUGS runs. 

When gas becomes cool (T < 15000 K) and 
dense {nth > 9.3 cm~'^), it is converted to stars 
according to a Kennincutt - Schmidt-like law with 
the star formation rate oc p^'^. Effective star for- 
mation rates are determined by the combination 
and interplay of and feedback, c^, is ultimately 
the free parameter that sets the balance of the 
baryon cycle off cooling gas, star formation, and 
gas heating. Stars feed energy back into surround- 
ing gas. Supernova feedback is implemented using 



the blastwave formalism (Stinson et al. |2006) and 
deposits lO'^^ erg of energy into the surrounding 
medium at the end of the stellar lifetime of ev- 
ery star more massive than 8 Mq . Energy feed- 
back from massive stars prior to their explosion 
as supernovae has also been included. Without 
the ability to resolve the details, we employ a rel- 
atively crude thermal implementation of radiation 
feedback from massive stars, with the aim being 
to mimic their most important effects on scales 
that we resolve, i.e. to regulate star formation, en- 
hance inhomogeneity, and to allow the expansion 
of the SNe driven super-bubbles which drive out- 
flows. To mimic the weak coupling of this energy 
to the surrounding gas ( Freyer et al.||2006" ), we in- 
ject pure thermal energy feedback, which is highly 
inefficient in these types of simulations (jKatzjl992] 
Kay et al.||2002|. We inject 10% of the available 



energy during this early stage of massive star evo- 
lution, but 90% is rapidly radiated away, making 
an effective coupling of the order of 1%. 

The two simulations have different initial mass 
functions, with G-1536 having a factor of ~ 2 more 



massive stars (Chabrier 2003) than G-1578411 
( Kroupa et al.|[l993 ), meaning that G-1536 is less 
efficient at turning baryons into stars. 

Ejected mass and metals are distributed to the 
nearest neighbor gas particles using the smoothing 



et al. 



kernel (Stinson et al. 2006). Literature yields for 
SNII ([ Woosley fc Weaver|19 95') and SNIa ( [Nomoto 
1997) are used. Metal are diff'used by 



treating unresolved turbulent mixing as a shear- 
dependent diffusion term (Shen et al. 2010), al- 



lowing proximate gas particles to mix their metals. 
Metal cooling is calculated based on the diffused 
metals. 

3. Bulge selection criteria 

The classical method of separating the contri- 
bution of the inner region from the disk rests on 
fitting the light radial profile of the spiral galaxy 
with two (or more) components (e.g. exponential 
disk -|- a Sersic profile). However, this method 
does not provide any way in which individual stel- 
lar particles can be assigned to the inferred galaxy 
components. Therefore we make a kinematic sep- 
aration of bulge and disk stars using clustering al- 
gorithms. For each star particle we computed Jp, 
the angular momentum in the plane of the disk, 
Jz, the angular momentum perpendicular to the 
disk, and Jc, the angular momentum of a parti- 
cle with the same binding energy (E), moving in 
a circular orbit at the same radius. Each stellar 
particle within the virial radius is then dynam- 
ically defined by three variables: Jz/Jc , Jp/Jc, 
and binding energy (E). These variables are nor- 
malized to [0,1] in order to give equal weight to 
all variables, and then fed to the clustering algo- 
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Fig. 1. — The Jz/Jc (left), Jp/Jc (center) and E (right) histograms for the galaxy G-1536. The spheroid 
and disk contributions, as derived with the clustering algorithm are given in red and blue, respectively. The 
histograms considering all the galaxy stellar particles are given in black. 




4 6 

R (kpc) 
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Fig. 2. — The projected stellar mass density for 
the kinematic disk (blue curve) and spheroid (red 
curve), and all stellar particles (black curve) of 
galaxy G-1536. The grey dashed line gives the 
position of the radial cut used to define the bulge 
region. 

rithm, which requires as prior only the number of 
clusters (we chose n^2 for bulge + disk). The 
normalization assumes a linear mapping between 
[Xjmnj'^max] and [0,1] (with J^/Jc, Jp/Jc and E 
as X) , which is appropriate if the histogram of X 
does not have extremely extended tails. We con- 
sidered the distance in the phase space to be the 
Euclidean one. The algorithm, starting from an 
initial random partitioning into n clusters, itera- 
tively searches for the partitioning which would 



minimize the intracluster distance. In Figure [T] 
we show the output of the clustering algorithm in 
the form of the histograms of Jz/Jc, Jp/Jc and E 
for disk, spheroid and galaxy stars in the case of 
G-1536. As expected, the spheroid stars peak at 
Jz/ Jc — in accordance with a system sustained 
by velocity dispersion, while the disk ones peak 
closer to Jz/Jc = 1- In the Jp/Jc histrograms the 
disk shows a peak at ~ 0.1, characteristic of a sys- 
tem whose particles move little perpendicular to 
the disk plane, although some of them can have 
an important vertical motion which shows up as 
the tail extending to Jp/Jc ~ 0.6. 

Once the dynamical decomposition has been 
done, the projected stellar mass density profiles of 
the disk (thin plus thick) and spheroid naturally 
provide us with the radial cut necessary to delimit 
the bulge, in the form of the radius, T:buige^ where 
the two intersect. Within this cut the mass contri- 
bution of the disk stars is much smaller than the 
spheroid, therefore minimizing the contamination 
from disk stars when determining the bulge global 
properties. An example of projected stellar mass 
density for the disk, spheroid and galaxy stellar 
particles is given in Figure [2] for G-1536, as well as 
the position of the radial cut we use to define the 
bulge region. 

Although a dynamical separation of galaxy 
components is of great help in analysis of galaxy 
structure formation, from an observational point 
of view the different dynamical components of 
the central galactic regions are usually hard to 
disentangle. For this reason, we also use an obser- 
vationally based selection for the bulge stars, in 
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the form of a simple radial cut at Thuige ■ 

Summing up, we use two bulge selection cri- 
teria: i) the stars within an sphere of radius 
fbuige, hereafter position selection, and ii) the stars 
within an sphere of radius rbuige that at the same 
time belong to the galaxy spheroid according to 
the kinematic-decomposition, hereafter kinematic 
selection. In Table [T] we give rj^uigei the total stel- 
lar and gas mass of each bulge at z = 0, as well as 
the mass resolution for each simulation. 

4. Bulge populations 

To analyze the bulge stellar populations and 
their link with dynamical processes, we con- 
structed Figure [3j in which both the star for- 
mation rate (SFR) histories of the bulge and the 
galaxy mass aggregation tracks (MATs) have been 
drawn. The MATs each give the evolution of mass 
inside a fix radius (right axis in each panel) . Virial 
radii have been calculated based upon the |Bryan| 
& Norman ( 1998 ) fitting function to determine 



the over-density threshold. For the baryonic com- 
ponent, radii are binned in equally spaced steps in 
a logarithmic scale. On the right (left) panels the 
colored curves show the stellar (stellar plus cold 
gas) masses within the respective radii as a func- 
tion of the Universe age {lu is the Universe age 
at z = 0) and redshift. MATs corresponding to 
radii within rij^igg^ are drawn with different colors, 
while those at larger radii are given in cyan. 

Major mergers {M secondary / M^ain > 0.25), mi- 
nor mergers and slow accretion processes in the 
dark matter or baryonic component can be clearly 
identified as big or small mass jumps in the MATs, 
or as continuous mass increments, respectively 
(note the different mass scales on the right axes for 
the P-DEVA and GASOLINE galaxies). Two differ- 
ent phases are reflected in the noticeable knee-like 
shape of the MATs in all objects shown in Fig- 
ure |3] an initial phase where the mass assembly 
rate is high, and then a much slower phase, when 
mass is more slowly acquired. 

Superimposed on the same figure are the bulge 
SFR histories (left axis in each panel) , for both po- 
sition and kinematic selected bulge populations. 
All objects present starbursts at large redshifts, 
peaking between z = 3 and z = \, where signifi- 
cant substructure is also apparent in many cases. 
These early star formation peaks are expected 



when the primordial gas suffers a violent collapse 
and are correlated with the knee-like feature vis- 
ible in their corresponding MATs. The high red- 
shift SF peak is a direct consequence of the first 
phase of mass assembly. The simulations also show 
later peaks in star formation that are related to 
merger/accretion events. 

At lower redshifts a variety of processes appear 
imprinted in the SFR histories: during quiescent 
periods, the SFR decays to an approximately con- 
stant tail at the same time as the MATs get flat- 
tened. The low tails present in all SFR histo- 
ries correspond to periods of mass assembly where 
only minor mergers or gas accretion at slow rates 
show up in the MATs. All simulation also show 
secondary jumps in SFRs, associated with later 
merger events. These vary in time, and in size 
relative to the initial starbursts, for example HD- 
5103B has a relatively late (t/t[/) accretion event 
that shows up in the total mass MAT (black line) 
as well as the build up of bulge mass and the star 
formation rate. 

Some qualitative differences are apparent in 
the two GASOLINE runs compared to the P-DEVA 
ones. The initial starbursts associated with the 
fast phase of mass accretion are less peaked in 
the former, due to the affects of feedback. In G- 
1578411, which has the lower mass but also lower 
feedback of the two GASOLINE simulations, the ini- 
tial starburst has a longer duration than the other 
simulations, while in G-1536 the initial starburst 
in the bulge region is more suppressed relative to 
the overall star formation within the bulge, due 
to the more efficient feedback in this simulation. 
Also, due to energetic feedback, the MATs of G- 
1578411 and G-1536 show delayed baryon mass as- 
sembly relative to their halos collapse at high zs. 
As a result, these objects have low central bary- 
onic mass density, the effect being more marked 
at high-z. 

Regardless of these specific differences, the cor- 
relations between MATs and SFR histories sug- 
gest that a meaningful stellar age classification can 
be based upon the two phases showing up in the 
MATs. Therefore we classify as old bulge stel- 
lar population the stars formed as a direct con- 
sequence of the fast phase (with formation times 
smaller than the temporal cut drawn in thick red 
in Figure [3| , while the stars formed later on we 
globally denote by young bulge stellar populations. 
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Fig. 3. — Bulge star formation rate comparisons for kinematic (solid gray) and position (solid black) selected 
bulge stars with sizes given by rijuige- The position of the thick red vertical (values in Table [2]) indicate the 
separation between the old and young stellar populations. The colored lines represent the mass aggregation 
tracks (MATs) of cold baryons (left) and stars (right) along the main branch of the merger tree for each 
object, the masses being computed inside fixed radii of 0.50, 0.72, 1.03, 1.47, 2.10, 4.30, 8.80 and 25.80 kpc 
from bottom to top. The upper solid black line is the total mass inside the virial radius (virial mass), the 
solid curves in red, green, blue, magenta and orange correspond to radii roughly within the bulge, while the 
dashed cyan ones correspond to radii exceeding r^uige (see Table [T]). 
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In a quantitative sense, the temporal cuts corre- 
spond to the time from which the second deriva- 
tive of the MAT^Ttuige) becomes fiat, reflecting a 
transition from a fast clumpy mass assembly to a 
slow smooth regime, without considering the MAT 
variations induced by low redshift mergers, like in 
the case of HD-5103B. The values of the temporal 
cuts, tcMt, for all bulges are given in Table [2] to- 
gether with the mass weight of the old population. 

5. 3D Shapes, sizes and kinematics 

To quantify the shape of the bulges, as well that 
of the old and young bulges as defined above, we 
show the correlations between the axis ratios of 
each ellipsoid of inertia in Figur e [4} The inertia 
tensor was computed following Gonzalez-Garci'a 



& van Albada ( 2005 ) , and subsequently diagonal- 
ized. Next, the eigenvalues (Ai > A2 > A3) were 
used to compute the length of the principal axis 
a ^ b > c. Red and blue represent the old and 
young bulge, while black gives the average over 
the whole object. Filled and empty symbols corre- 
spond to the objects in the position and kinematic 
selection, respectively. 
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Fig. 4. — The axis ratios of the five spheroids. 
Red and blue represent the old and young bulge, 
while black gives the average over the whole ob- 
ject. Filled and empty symbols correspond to the 
objects in the position and kinematic selection, re- 
spectively. The dashed-dotted gray curves sepa- 
rate prolate, triaxial and oblate morphologies. 
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Fig. 5. — Minor vs major semi axis for the five 
spheroids. Red and blue represent the old and 
young bulge, while black gives the average over 
the whole object. Filled and empty symbols cor- 
respond to the objects in the position and kine- 
matic selection, respectively. The dashed line is 
where c=a. Young bulge populations are smaller 
than old ones, with the exception of G-1578411. 

in Figure |4] have h — c, while the dashed-dotted 
gray curves denoted by T = 0.3 and T = 0.7, 
respectively, separate the oblate objects from the 
triaxial and the prolate ones according to the T pa- 



rameter introduced by de Zeeuw & Franx (1991) 
with the definition T = (1 - (&/a)^)/(l - [cjaf)- 
The oblate objects correspond to c/a < 0.9 and 
T < 0.3, the prolate to c/a < 0.9 and T > 0.7, 
and the rest to triaxial ones. In this perspective, 
we observe that both complete bulges as well as 
their distinct components are all oblate with the 
exception of the secular bulge of G-1578411 which 
is triaxial regardless of the selection criteria. Old 
and young populations are clearly segregated, ei- 
ther in the position or in the kinematic selections, 
and more separated in c/a than in hja. Old popu- 
lations are more spheroid-like, while the young are 
more oblate, consistent with observations. 

Figure [5] gives an idea of the sizes of the five 
bulges by plotting c vs a. In all but one case (the 
a value of G-1578411), the old bulges have larger 
sizes (as given by a and even more by c) than 
the young ones, sustaining the idea that younger 
stellar populations ocupy a smaller volume than 
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Table 2: The temporal separation between the distinct star formation episodes for the simulated bulges, the 
corresponding mass percentage of old and young stars in the position-selection, and the Sersic indices 
derived from fitting the projected stellar mass density of old and young (according to the bulge temporal 
cut), and of all galaxy stellar particles. 




Fig. 6. — Rotational support versus ellipticity. 
Red and blue represent the old and young bulge, 
while black gives the average over the whole ob- 
ject. Filled and empty symbols correspond to the 
objects in the position and kinematic selection, 
respectively. Young bulges are more rotationally 
supported and less spherical than the old popula- 
tions, even when kinematically selected to exclude 
the disk. 

their corresponding old counterparts. In this re- 
spect the following inequality holds, irrespective of 
the selection criteria used: y/ UybyCy < \/atbtCt < 
\/aoboCo, where y, t and o stand for young, total 
and old, respectively. This graph shows a clear 
sequence, the size increasing from the younger to 
the older, with the total population in between. 
Also, as it was proven by the previous graph, the 
older bulge components are closer to be perfect 
spheroids (they appear closer to the dashed line 
for which c — a). 



The difference between the results obtained us- 
ing the two selection criteria are only qualitative. 
Indeed, the trends in shapes and sizes are the same 
for both selections. The kinematic-selection leads, 
on average, to bulge components closer to being 
spherically symmetric than the position selection 
does. This is an expected result since position 
selection implies the exclusion of precisely some 
of the particles encountered in ordered rotation 
which would otherwise lead to ellipsoids of inertia 
characteristic of more disk-like objects. 

Also, the shape and size trends with stellar age 
are the same for the P-DEVA and GASOLINE bulges. 
The only difference is that GASOLINE bulges have 
lower mass concentrations and therefore larger val- 
ues for a and c, due to energetic feedback. 

Let us now turn to the bulges kinematical anal- 
ysis, including some of their classical or pseudo- 
bulge properties, therefore using line-of-sight (los) 
velocities. In order to analyze the rotational sup- 
port, we align the galaxies with the z-axis per- 
pendicular to the disk. Since the rotational ve- 
locity can be best observationally measured in 
edge-on projection, we consider two lines of sight 
(along the x- and y-axes), and compute the ra- 
dial profiles of the velocities - wios=x(^ = y,z) 
and ■i;ios=y(^ = x, z), respectively - along the two. 
The dependence of the profiles on the altitude z 
above the galactic plane is important both when 
studying whether the bulge or its components are 
in cylindrical rotation, as well as from the point 
of view of comparison with observations. In the 
latter case, data are normally taken, due to the 
extinction effect in the disk plane, with long-slit 
or IFUs at altitudes higher than the vertical scale 
height of the thin disk. Considering all these ob- 
servational limitations, we constructed the veloc- 
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ity profiles using a 2;-binning of 0.2 kpc, averag- 
ing the curves in the two hemispheres in order 
to minimize the statistical noise in the slits at 
higher latitudes where the number of particles is 
smaller. Complementary to the profiles of rota- 
tional velocity, we also constructed the los veloc- 
ity dispersion radial profiles, aios{R,z), which are 
approximately flat. Given the flatness of these 
profiles, we considered the crios(2) as the average 
of crias=xwiR — y I 2;, z) over all lines of sight 



R = y 



and weighted with the number of par- 



ticles in each bin nu- Once the rotational velocity 
profiles were constructed, we took as Vmaxiz) the 
weighted average of the profile extremes. Finally, 
we averaged Vmax{z) and aios{z) over the two lines 
of sight, los=a; and los^y. 

In order to choose a representative z-bin for 
all bulges, we considered both the disk vertical 
scale heights as well as the bulge radial extensions. 
Therefore, we selected the bin with 0.4 <| z |< 0.6 
kpc. Given this choice, within the current section, 
Vmax and (Tios will refer to this specific slit. 

Figure [6] plots rotational support, defined as 
the maximal rotational velocity (vmax) divided by 
the bulge los velocity dispersion (crios), as a func- 
tion of the intrinsic ellipticity e — 1 — c/a. The 
dashed black curve in the graph represents oblate- 
spheroid systems with isotropic velocity disper- 
sion flattened only by rotation ( Binney|1978| ) , and 
is approximately described by y'e/(l — e) (Kor- 
mendy||1982] ). 



The results show a clear trend of higher eccen- 
tricities and rotational support for the younger 
populations (blue symbols) as compared to the 
older ones (red symbols), while the whole bulges, 
in black, occupy intermediate positions. It is im- 
portant to note that these patterns hold irrespec- 
tive of the code used to run the simulations or the 
bulge selection criteria, as well as of the slit po- 
sitioning. Note however that increasing (decreas- 
ing) the z slit positioning lowers (increases) the 
specific values of Vmax/crios for all the five bulges 
as well as in their distinct stellar components. 
Therefore, our bulges appear more classical as we 
move away from the galactic plane, but always 
they show same trends with stellar age, with the 
younger component being more rotationally sup- 
ported. We note that in all cases, a weighting 



by light rather than mass will provide increased 
prominence to the younger pseudo-bulge-like pop- 
ulations. In two cases (HD-5004A and HD-5103B 
in the dynamical selection), the bulge as a whole 
has a slightly lower value of Vmax I <^\os even than 
the old population. In these cases, the bulges have 
velocity dispersions approximately equal to the old 
components, but have slightly smaller Vmax- This 
difference in v^^ax can be explained through mis- 
alignments and/or kinematical peculiarities. 

Another parameter describing the bulge shape 
is the Sersic index [n). Thus, we fitted the pro- 
jected stellar mass density profiles of our five 
galaxies with a Sersic -I- an exponential disk. The 
data and the fits are depicted in Figure [7| while 
the Sersic indices are given in Table [2] With the 
aim of checking whether n varies if considering 
only the old or only the young stars, we also fit- 
ted the two stellar subpopulation of the galaxies, 
using for the temporal cut the limit defined with 
respect to the bulge region. It is important to 
stress, that n values are to a large extent depen- 
dent on the amount of feedback. In the case of 
GASOLINE galaxies, the incresed feedback leads to 
a fiattening of the projected stellar mass curves at 
small radii. For this reason, the Sersic indices of 
G-1578411 an G-1536 are < 2 irrespective of con- 
sidering all galaxy stars or only the corresponding 
components, while those of HD-5004A, HD-5004B 
and HD-5103B (simulated without explicit ener- 
getic feedback) are > 2. In any case, the differ- 
ence between rioid and Uyoung is small in any of the 
simulated galaxies, with a slight tendency for the 
old stars to have larger n than the younger ones, 
with the exception of G-1536, which was simulated 
with a considerable amount of feedback. From an 
observational point of view, no define trend with 
the band has been detected either, see for example 



Fisher & Drory (2008), or more recently McDon- 



ald et al. (20111. The Sersic index has been used. 



together with other parameters like rotational sup- 
port, to classify bulges. For example, |Kormendy fc| 



\ X stands for y or x 



Kennicutt (2004) consider pseudo bulges to have 
< 2, while classical ones have > 2. In this respect, 
the three P-DEVA galaxies have classical bulges, 
the old bulges being more classical than the young 
ones {rioid > Uyoung)- On the other hand, the two 
GASOLINE ones have pseudo bulges according to 
this classification, and show no definite trend with 
age. In the same perspective, the pseudo bulges of 
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Fig. 7. — The Sersic + exponential disk fits of the projected mass density of galaxy stars (left panels), and of 
the corresponding old (central panels) and young (right panels) stellar sub-populations, separated according 
to the bulge temporal cut in Table [2j The color code is: grey open squares for the data, red and blue lines 
for the Sersic and exponential disk contributions, and black lines for the total fit. 
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Okamoto (20121, formed at high redshift, roughly 



during what we define as the fast phase of mass 
assembly, have < 2. However, his simulations also 
use energetic feedback which is at least partially 
responsible for the low values of n. Therefore, we 
think that further studies are needed in order to 
disentangle to what extent the n values depend 
on the particular modes of mass assembly and to 
what extent on the amount and implementation 
of the feedback effects. In this respect, our study 
suggests that, at least for simulations, v/a sepa- 
rates more reliably the young from the old bulge 
than n, which tends to be similar for old and young 
stellar populations. 

6. Ages and abundances 

This section concerns the ages and metallicities 
of the old and young bulge stars at z = 0. 

Although details differ slightly, in both the 
P-DEVA and GASOLINE simulations trace the pro- 
duction and enrichment of chemical elements in 
broadly the same way. Newly produced elements, 
as by-products of stellar evolution and death, are 
released to the surrounding interstellar medium as 
increments in the metal content of nearby gas par- 
ticles. Metal diffusion is implemented among gas 
particles with a diffusion constant, allowing for el- 
ements to mix in a given environment, mimicking 
what happens in a turbulent interstellar medium. 
Both codes consider the evolution of the follow- 
ing elements: H, He, C, N, O, Ne, Mg, Si, S, Ca 
and Fe. As we are looking at broad trends here, 
the relatively small differences in the implementa- 
tion chemistry between the codes are not critical. 
As a measure of metallicity we used either [Fe/H] 
or the average values of these elements; also, we 
used [Mg/Fe] to follow a-element properties and 
checked that using [0/Fe] does not change sub- 
stantially our results. The reference values for 
the solar metallicities were taken from IGrevessel 



& Sauval (1998). 



We analyze the distributions of metallicity and 
a-element abundances in the different bulge pop- 
ulations. They are plotted in Figure [s] (central 
and right panels), where it can be seen that the 
abundance distributions of young and old bulge 
populations are segregated (the old stars have a 
lower metallicity and a higher [Mg/Fe]), giving in 
most cases overall two-peak [Fe/H] and [Mg/Fe] 



distributions. The slightly cleaner separation of 
two populations in [Mg/Fe] space in the P-DEVA 
simulations as compared to the GASOLINE ones 
might partially be an effect of the energetic feed- 
back in the latter, which delays the SF relatively 
to dynamical processes (see Figure [s] and the cor- 
responding comments in Section |4| and enhances 
the mixing, although the differences in the de- 
tails of implementation of metal enrichment could 
also play a role, as well as in the wider ranges 
of [Mg/Fe] of the P-DEVA runs. However, the 
important point here is that both codes result in 
broadly the same relative trends in metal abun- 
dances and enrichments of their old and young 
populations. 

In this respect, we also plot [Mg/Fe] as a func- 
tion of [Fe/H] (left panel). Again the young and 
old populations show up at different loci in this 
plot. A feature to be noted in the graphs [Mg/Fe] 
versus [Fe/H] are the slopes of the two stellar pop- 
ulation types in the sense that the old one has a 
milder slope, while the young population shows a 
steeper one. Again, these trends are apparent in 
all simulations, and what is interesting is that the 
age selection we have used invariably separate at 
the "knees" of the abundance trends. This feature 
has to do with the different nucleosynthetic origin 
of a and Fe elements, a clean knee- like pattern be- 
ing expected in stellar populations where the SFR 
is concentrated in an episode with a very short 
time-scale followed by an epoch of more quiet star 
formation. The plateau corresponds to the early 
stages when SNe II dominate the metal enrich- 
ment, followed by a downturn to lower [Mg/Fe] 
when enough time has elapsed since the high red- 
shift starburst that SNe la come into play (see 
Wyse|[l999 and references therein). 



These results on the chemical composition of 
bulge populations agree with recent observations 
of the Milky Way Bulge, like those of |Babusiaux| 
et al. (2010), who found two distinct populations 



in Baade's Window, one old, metal-poor compo- 
nent with kinematics of a spheroid and a younger, 
metal-rich one displaying kinematic characteristics 
similar to a bar. Otherwise, Ness et al. ( 2012[ ) ex- 
tended this work by enlarging the number of win- 
dows, finding consistent results (see Section [l] for 
more details). 

Figure |9] gives the averaged metallicities and 
ages of the distinct stellar populations as functions 
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Fig. 8. — [Mg/Fe] vs [Fe/H] (left), and the [Fe/H] (center) and [Mg/Fe] (right) histograms lor the simulated 
bulges. The bin widths for [Fc/H] and [Mg/Fe] are 0.07 and 0.04 dex, respectively. Red and blue represent 
the old and young bulges. 
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Fig. 9. — Average metallicity (top) and age (bot- 
tom) of the distinct bulge components versus los 
velocity dispersion. Red and blue represent the 
old and young bulge, while black gives the aver- 
age over the whole object. Filled and empty sym- 
bols correspond to the objects in the position and 
kinematic selection, respectively. 



of the los velocity dispersion. In this case, crios was 
computed in the same way as in SectionjS] the only 
difference being that no z-binning was used. First 
of all, we note that a clear population segrega- 
tion stands out in both panels, irrespective of the 
code or the selection criteria employed. In the up- 
per panel, the metallicity appears to be increasing 
with the velocity dispersion, when considering all 
five bulges. However, if discarding the G-1578411 
bulge (bottom- up triangles), {Z) seems to be al- 
most constant with crios. On the other hand, if 
considering only the P-DEVA bulges, (Z) decreases 
with (Tios. These trends should be considered with 
caution, given the inhomogeneity of our sample. 
In the age - cios plot, we get a positive slope, in 



accordance with the findings of MacArthur et al.] 
(|2009|). Taking dios as a measure of the objects 



mass ( Binney & Tremaine 1987 ) (see Table [I] for 
the bulge stellar masses), it can be noted that the 
more massive objects also have the older average 
overall populations. 

7. Bulges within The Cosmic- Web 

We trace the evolution of particles in the bulge 
at low ziow = to their origins at Zj„ = 10, fol- 
lowing them at 60 different redshifts in between. 



In Figure 10 we plot the evolution of the baryons 
that form the bulge stars of IID-5004A at z = 
(left and central graphs) at four zs (top to bot- 
tom) representing relevant consecutive events in 
this bulge formation, as well as the object's sur- 
roundings (right graphs) at each redshifi]^ In 
the other simulations this sequence of evolution 
is similar enough that we can use this as a rep- 
resentative example. The left sequence aims at 
illustrating the differences in the space configura- 
tions that old z = bulge-to-be stars (red points) 
form along their evolution as compared to those 
of young z — bulge-to-be stars (blue). The 
emphasis here is on the two population types, ir- 
respective whether they are gas or stars at the 
Universe age corresponding to the plot. The cen- 



^Indeed as Figure [s] shows, the first snapshot at t/tu = 
0.13 gives us the configuration of the bulge-to-be stars at 
the very beginning of the collapse-like event; the second 
one at t/tu = 0.28 roughly corresponds at the end of the 
fast phase; the third one at t/tu = 0.38 is just before the 
beginning of the major merger causing the strong stellar 
burst in the slow phase, and finally the fourth one a,t t/tu = 
0.59 corresponds to the end of this burst. 
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tral sequence takes into consideration the bulge- 
to-be stellar particle properties at the zs plotted, 
namely points are green for cold gaseous particles 
(T(z) < lO'^K), black for hot ones (T(z) > 10^ if), 
and red and blue for stellar particles, depending 
on their time of birth according to Table [5] Its aim 
is to illustrate the differences in the birthplaces of 
the two bulge stellar populations relative either to 
the cosmic-web (at high zs) or to the protogalaxy 
components (at lower zs) , as well as their evolution 
towards their final configuration at z = 0. Finally, 
in the right panels we plot all the particles in the 
same volumes, and with the same color code as in 
the central one, in order to compare the dynamical 
processes of the bulge-to-be stellar particles with 
those of the surroundings mass elements causing 
them. 

We first analyze the sequence of bulge forma- 
tion events in terms of the cosmic-web element 
dynamics (central panels). Before z = 3.5 mass 
piles up in caustics |^ where already some nodes 
show up. At this redshift, star formation is al- 
ready triggered in the densest of these nodes at 
separated places. Comparing the two upper snap- 
shots in any of the vertical sequences, we see that 
an overall collapse-like event acts in between onto 
a structured net of cells as a contractive deforma- 
tion. It somehow erases the cell structure, join- 
ing together clumps mainly along filaments, and 
therefore only very low relative angular momen- 
tum is involved at these zs. The collapse-like 
event shrinks the volume visualized at z = 3.5 into 
that at z = 1.75, where a central mass concentra- 
tion fed by filaments stands out. Gas keeps on 
flowing through filaments towards the central re- 
gions between snapshots two and three. A diffuse 
component (i.e., outside caustics) is also evident 
at high z. This diffuse component tends to flow 
into caustics, therefore vanishing as evolution pro- 
ceeds. Otherwise, by z = 1.21 the filaments have 
practically been removed in favor of clumps and 
some cold gas in irregular structures. Moreover, 
two small gaseous disk-like structures and a third 
central disk with a stellar component, seen edge- 



^ Caustics are the elements of the cellular structure where 
dense mass elements show up at high zs. They are classified 
into walls, filaments and nodes, and, at a given scale, these 
represent a temporal sequence of mass piling up. See|Shan- 
|darin Sc Zeldovich] (jl989 ) and |Dommguez-Tenreiro et al, 
^2011^ for" more details. 



on, have formed in the interval between the second 
and third snapshots. At z = 1.21, the stars are at 
the center of these disks, or at the center of other 
smaller gaseous structures. The young-to-be bulge 
stars begin to show up (blue points). Finally, after 
a major merger occurring between the third and 
fourth snapshots, by z = 0.58 practically all the 
baryonic particles that are to form the HD-5004A 
bulge stars at z = come to be bound in a unique 
system (except for some small clumps). Most of 
the particles that have transformed into stars (red 
and blue at the central sequence) are at the center, 
while those that are still gaseous (in green) show 
up in a rebuilt disk-like configuration. This config- 
uration disappears later on, and by z = 0.35 (not 
shown in Figure 10 ) practically all the bulge par- 



ticles are at their place within a sphere of radius 

^bulge • 

It is very useful to compare these evolutionary 
processes with the mass distribution of their sur- 
roundings visualized in the right graphs of Fig- 
This comparison makes it even clearer 
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that bulge material has been involved in caustic 
formation at high z although part of it remains 
as diffuse gas quite a while. Interestingly, a frac- 
tion of it has formed transient gaseous disks or has 
been a part of the former disk of the HD-5004A 
object. 

We now analyze the differences in the assembly 
patterns of the young and old bulge stellar popu- 
lations relative to the cosmic-web dynamics (left 
graphs in Figure 10). Most mass elements identi- 
fied at z = as old stellar particles (red points) 
are already concentrated in caustics by z — 3.5. 
These are the densest regions at this time. In con- 
trast, an important fraction of the mass elements 
identified at z = as young stellar particles (blue 
points) are not yet involved in caustics. The differ- 
ence increases as evolution proceeds. In fact, the 
strong contractive deformation acting between the 
first and second snapshots involves preferentially 
the old-stellar-to-be particles, in such a way that 
by z = 1.75 the old component has a very small 
volume, while the young component still spans 
a scale of ^ 200 kpc and has an important dif- 
fuse gas fraction (with some nodes formed at lo- 
cal contractions as well). In the third snapshot 
we see that a fraction of the mass elements that 
are to form young stars now shows disk-like pat- 
terns, while old stars are at the centers of these 
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Fig. 10. — The different components of the bulge HD-5004A (left and central panels) and its surroundings 
(right panels) at redshifts 3.50, 1.75, 1.21 & 0.58 {t/tu of 0.13, 0.28, 0.38 & 0.59), from top to bottom. The 
left panels give the positions of the bulge stellar particles identified at z = when traced back to each of 
the four zs. Red corresponds to baryon particles (either gas or stars) whose transformation into stars occurs 
along the fast phase of bulge formation (i.e., at tform/if/ < 0.31, or equivalently, z > 1.54 in this particular 
case), and blue to bulge baryon particles that become stars at later times. The same traced back positions 
are given in the central panels, but in this case with a color code representing particles properties at the zs 
plotted, namely green for cold gaseous particles (T < lO^K), black for the hot ones, and red and blue for the 
stellar particles according to their time of birth (older and respectively younger than z = 1.54). The right 
panels show the positions of all baryonic particles within a limiting radius of the galaxy's progenitor center 
using the same colors as in the central graphs. 
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small disks. These disk-like configurations pro- 
vide a key to explain the different kinematic prop- 
erties of young bulge populations as compared to 
old ones, because of the angular momentum con- 
tent of the former. As said above, the disks seen in 
the third snapshots merge. This causes their cen- 
tral old stars to form a unique bulge- like structure. 
At the same time, new stars - belonging to the 
young bulge - form and appear within it. On its 
turn, the young component, partially still as dif- 
fuse gas, forms structures (i.e. a disk in the case of 



HD-5004A as it can be seen in Figure 10 ) around 
the old stars spheroid. Some of the small nodes 
remain as satellites, carrying a small bit of old 
population at their centers. Finally, by z = 0.35 



(not shown in Figure 10), both the old and young 
bulge populations (or the gaseous particles to form 
the latter) are bound to their respective central 
spheroids, where they remain until z = 0. 

We draw attention to similarities with mas- 



sive ellipticals Dommguez-Tenreiro et al. (20111 



In this case, stellar particles at z = exhibit a 
gaseous web-like morphology at z ~ 3.5 6, with 
scales of ~ 1 physical Mpc. The densest mass ele- 
ments of this gaseous web, dynamically organized 
as attraction basins for mass flows, have already 
turned into stars by z ~ 6. At high z these basins 
undergo fast contractive deformations which vi- 
olently shrinks them in quasi-radial directions, 
therefore involving very low angular momentum. 
They can be described as collapse events with very 
complex geometries, causing high rates of dissipa- 
tion and stellar formation out of the available gas, 
most of if it being transformed into stars at the 
end of this process. Afterwards, during the second 
phase, the mass assembly rate is much lower and 
is characterized by mergers involving significantly 
larger amounts of angular momentum. 

To sum up, the old and young populations of 
the bulge HD-5004A are not only segregated in 
their z=0 properties, but also show quite different 
assembly patterns in terms of the cosmic-web dy- 
namics. Overall, this scenario also holds for the 
other four bulges. Indeed, most of the old stars 
formed at disjoint places within attraction basins 
during the early fast processes, and undergo an 
important assembly episode later on, through mul- 
ticlump collapse-like events shrinking these basins 
in quasi-radial directions. This assembly phase re- 
sults in a relative low angular momentum of the 



merging clumps. 

The ancestors of young bulge stars, on their 
turn, in many cases have been part of disk struc- 
tures before arriving to the galaxy central regions. 
Different destabilizing processes can be respon- 
sible for this inward material transport, such as 
mergers or secular processes in the disks. Young 
bulge stars keep partial dynamical memory of their 
angular momentum content. These differences in 
the assembly patterns relative to the angular mo- 
mentum involved, would be a key piece in ex- 
plaining the segregated properties of the two bulge 
components. Although our sample is small, the 
mass weight of the old component varies between 
31 and 50%. An even more important variation 
of this weight, together with the different inward 
mass transport mechanisms in the second assem- 
bly phase, could explain the important dispersion 
in bulge properties. 

8. Summary, Discussion and Conclusions 

We have analyzed the bulges of an inhomoge- 
neous suite of five spiral galaxies emerging from 
high resolution hydrodynamical simulations in a 
cosmological context. Our aim is to decipher the 
underlying physical processes in bulge assembly 
that could be responsible for their observed prop- 
erties, focusing on relevant structural, kinemati- 
cal and chemical properties of their stellar pop- 
ulations and paying particular attention to their 
observationally suggested segregation by age. 

Concerning the problem of bulge formation, the 
main particularity of this work is that we ana- 
lyze simulated disk galaxies that are run with dif- 
ferent codes (including different sub-grid physical 
prescriptions) where the chemical evolution has 
been carefully implemented, by using a common 
pipeline to measure the relevant stellar popula- 
tion properties. The simulations have been run 
with P-DEVA and GASOLINE. Feedback from mas- 
sive stars is implemented explicitly in the two 
GASOLINE simulations through the blast-wave for- 
mulation (IStinson et al.l 120061), while in P-DEVA 



simulations the effects of discrete energy injection 
are assumed to be on subgrid scales, resulting in 
the low star formation efficiency which is assumed 



to mimic them (see discussion in Agertz et al. 
20TT| ). 

Although details differ slightly, both P-DEVA 



19 



and GASOLINE simulations trace the production 
and enrichment of chemical elements in broadly 
the same way. The newly formed elements, as 
by-products of stellar evolution and death, are re- 
leased to the surrounding interstellar medium as 
increments in the metal content of nearby gas par- 
ticles. The metals are diffused among gas parti- 
cles, allowing the elements to get mixed and there- 
fore, mimicking the turbulent interstellar medium. 
Both codes consider the evolution of 11 elements. 
As we are looking for broad trends here, the rel- 
atively small differences in the implementation 
chemistry between the codes is not critical. 

In this paper we have studied in detail the mass- 
averaged three-dimensional sizes, shapes and kine- 
matics, as well as stellar ages and element compo- 
sitions of these five bulges. This is an intrinsic ap- 
proach as opposed to looking for quantities closer 
to observations (i.e., light- averaged). We keep as 
this level because our aim in this paper is to find 
out the patterns of bulge mass assembly by focus- 
ing on their stellar population properties. 

We have found a satisfactory qualitative agree- 
ment with the latest observational data for the 
Milky Way as well as for other external bulges 
(see Sectionjl]). Our results indicate that bulges in 
our sample have an old stellar population, formed 
at high z at disjoint places within attraction 
basins, and joined together through their rapid 
quasi-radial multiclump collapse, where the rela- 
tive angular momentum of the collapsing/ merging 
clumps is low. A second phase follows, with lower 
mass assembly and star formation rates, but with 
higher angular momenta. This phase shows a va- 
riety of assembly patterns: i) only minor mergers 
without further significant SF bursts but a SF tail, 
ii) major mergers with secondary SF bursts, as 
the galaxy HD-5004A shows in Figures [3] and 10 



and/or iii) secular evolution of the galactic disk, 
as is the case for G-1578411, which was shown in 
an earlier study to form its young bulge through 
secular processes after z = 1, at least partially 
driven by a bar ( Brook et al.||2012a ). 

The sizes, shapes, kinematics, stellar ages and 
metal contents of the stellar populations formed in 
these distinct phases, can be nicely distinguished 
in simulated bulges. Indeed, we have found that 
the young component tends to occupy a smaller 
volume, to have disk-like morphology (note that 
G-1578411 is rather triaxial), to be more rotation- 



ally supported, to have roughly solar metallicities 
and sub-solar a-element enhancements. The old 
population, by contrast, is more spheroid- like, has 
sub-solar metallicities and larger a-element en- 
hancements. On the other hand, no clear trend 
with age shows up in the Sersic indices, a re- 
sult in agreement with observations. The stel- 
lar metal content as well as the [Mg/Fe] ratios of 
these bulges have segregated distributions and, in 
some cases, show two clearly distinguishably peaks 
corresponding to the old and young populations. 
These stellar populations are also clearly segre- 
gated by their loci in the [Mg/Fe] vs [Fe/H] plots, 
where the old population has only a mild slope, 
while the slope of the young population appears 
steeper, with the two populations meeting at the 
knee. This kind of behavior is expected given the 
different nucleosynthetic origin of a and Fe-group 
elements. 

These trends have been shown to be robust 
against the different codes used, which have dif- 
ferences in their gravitational and hydrodynamical 
force integrators, in mass and spatial resolutions, 
simulation box sizes, star formation parameteriza- 
tions, chemical feedback and evolution, and ener- 
getic feedback implementations. Also, by chang- 
ing the bulge identification criteria from a simple 
radial cut to a kinematic based one, did not affect 
the age tendencies described above. The temporal 
cut-offs used to separate the two stellar popula- 
tions have been chosen based on the shapes of the 
galaxy mass aggregation tracks and on the behav- 
ior of the SFRs at the bulge scale. Varying slightly 
these cut-offs does not affect these tendencies, ei- 
ther. 

If we associate the old populations formed dur- 
ing the rapid phase with classical bulges, and the 
young ones formed during the slow phase with 
pseudo bulges, all the simulations in this paper 
show both classical and pseudo bulge components, 
but with varying relative masses. Therefore, the 
classical vs pseudo-bulge characteristics would be 
a question of degree, rather than nature. Other- 
wise, in all cases, measuring the light will result 
in higher relative contributions of the slow phase 
of bulge formation, i.e., pseudo bulges. In this 
respect, we note that the feedback in G-1536 is 
particular effective at quenching star formation in 
the early phase, resulting in a relatively less sig- 
nificant (in mass) classical bulge population (see 
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Table This may be pointing toward the solu- 



tion to the issue raised in Kormendy et al. (20101, 



who point out that a significant number of local, 
massive spiral galaxies have pseudo bulges rather 
than classical bulges. 

These assembly patterns are reminiscent of the 
two phases found in hydrodynamical simulations 



by Dommgucz-Tenreiro et al. (20061; Oser et al 



( 2010) ; Domhiguez-Tenreiro et al. (2011) for more 
massive early type galaxies. The main difference 
lies in the percentage of gas transformed into stars 
at early epochs. In the case of massive ellipticals 
most of the available gas at the attraction basins 
for mass flows is transformed into stars during its 
"collapse" along the first phase. On the contrary, 
no such exhaustive gas consumption occurs for the 
less massive galaxies, where gas remains available 
along the slow phase, first as a diffuse component, 
and later on being part of different structures (the 
host galaxy disk itself, other small disk-like satel- 
lites and/or small clumps) before it is incorporated 
to the bulge. 

If indeed real bulges follow a similar forma- 
tion pattern as simulated ones, their observational 
features discussed in Section [T] can be nicely ex- 
plained. More so, this approach provides a pos- 
sible explanation for some apparently paradoxical 
observational results. For example, metal gradi- 
ents could result from the different space distribu- 
tions of the old and young populations, with the 
latter being more concentrated at the center due 
to dissipation. Also, bulge rejuvenation can be 
easily explained within this scenario. 

We conclude that bulges can follow different 
assembly patterns, which can be summarized as 
two-phase processes (as in ellipticals) where non 
exhaustive gas transformation into stars occurs in 
the fast phase (unlike in ellipticals), with the addi- 
tional effects, along the slow phase, of major and 
minor mergers, as well as of disk secular instabil- 
ities, in some cases. These different patterns and 
their combinations in different epochs as well as 
in different proportions, might explain the impor- 
tant dispersion in bulge properties observationally 
found. 

The assembly of bulges is driven to a large ex- 
tent by dynamical processes at larger scales. This 
is particularly true concerning bulge-forming star- 
bursts. As our simulations show, the bulge mass 
aggregation is a delayed consequence of its host 



galaxy mergers in the slow phase, and the result of 
collapse-like events in the fast one. We have shown 
that galaxy mass and feedback can affect the rela- 
tive contributions of these two phases, although it 
is certain that specific basin deformation/collapse 
processes, as well as merger histories will also play 
a role. A set of statistical samples of simulation 
runs with various physical parameters is needed 
in order to break the degeneracy with the merg- 
ing history and provide further insights into the 
importance of the two phases of bulge formation. 
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